function sb = solve_backward_counter(param, fp,fp_parfor)

alpha = param(1:9);
TFP = param(10:11);
delta = param(12:14);
gamma = param(15:19);

%For Parents who are moved into the new neighborhood (they keep their preferences);
param_preference_PS_moved= delta(2) + param(20)*(fp.mean_income_by_neighborhoods(fp.origin_neigh) - fp.mean_income_by_neighborhoods(fp.receiving_neigh) );
%Discount Factor;
beta = 0.95;

estim_policy0=cell(fp.periods-1,2);
estim_policy1=cell(fp.periods-1,2);

estim_value  =cell(fp.periods-1,2);
estim_value0 =cell(fp.periods-1,2);
estim_value1 =cell(fp.periods-1,2);

%Grid points;
grid_peers=fp_parfor.grid_peers;
inv0 = fp_parfor.inv0;
inv1 = fp_parfor.inv1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Create Grid;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

hk = fp_parfor.grid ;

if fp_parfor.i==1
    %Initial guess in the solution algorithm;
    Hc_tp1_guess{1} = max(fp_parfor.distrib_skills(:,2)).*ones( length(inv1{1}), length(hk{1}) , length(grid_peers{1}) , 2 );
    Hc_tp1_guess{2} = max(fp_parfor.distrib_skills(:,3)).*ones( length(inv1{2}), length(hk{2}) , length(grid_peers{2}) , 2  );

else

    % Previous distribution in the fixed-point algorithm;
    Hc_tp1_guess = fp_parfor.Hc_tp1_guess_prime ;

end



for t=fp.periods-1:-1:1

    grid_p=grid_peers{t};
    grid_p = sort(grid_p);


    if t==fp.periods-1     %Solve last period

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Solving the parent's problem;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        grid_inv0 = inv0{t};
        grid_inv1 = inv1{t};

        %P_types stands for whether parents is moved or not. Moved parents
        %keep their preferences for parenting;
        [I0, hc0, H0, P_types0] =ndgrid(grid_inv0, hk{t}, grid_p, [0,1] );
        [I1, hc1, H1, P_types1] =ndgrid(grid_inv1, hk{t}, grid_p, [0,1] );

        log_Skills_tp1_0 = log(technology(t,alpha, hc0, I0, 0, H0, TFP) ) ;
        log_Skills_tp1_1 = log(technology(t,alpha, hc1, I1, 1, H1, TFP) ) ;
        
        %Parent's problem in the last period (equation 3) ;
        U0 = delta(1).*log(hc0).*(1 - delta(3).*0) + 1.*log(1-I0)  -(param_preference_PS_moved.*(P_types0==1)+delta(2).*(P_types0==0)).*0  + beta.*delta(1).*log_Skills_tp1_0;
        U1 = delta(1).*log(hc1).*(1 - delta(3).*1) + 1.*log(1-I1)  -(param_preference_PS_moved.*(P_types1==1)+delta(2).*(P_types1==0)).*1  + beta.*delta(1).*log_Skills_tp1_1;

        [Value0,ind0] = max(U0,[],1);
        [Value1,ind1] = max(U1,[],1);


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Computing Value e Policy Functions;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;


        for ps = 1:1:2 %Loop over types of parents (moved or receiving);

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
            %Value Function;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
            Value0_actual(:,:,:,ps,t) = Value0(:,:,:,ps);
            Value1_actual(:,:,:,ps,t) = Value1(:,:,:,ps);


            %Collect State Variables;
            hc_s = hc0(1,:,:, ps);
            H_s  = H1(1,:,:, ps);
            State_var = [ hc_s(:), H_s(:) ];

            %Value for P=0;
            value0 = Value0(:,:,:,ps);           
            Y0 = value0(:);
            estim_value0{t, ps} = compute_value_function(Y0, State_var);    

            %Value for P=1;
            value1 = Value1(:,:,:,ps);
            Y1 = value1(:);
            estim_value1{t, ps} = compute_value_function(Y1, State_var);    

             %Expected Value over preference shocks;
             Y = log( exp(Y0) + exp(Y1));
            estim_value{t,ps}  = compute_value_function(Y, State_var);


            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Policy Function (Investments)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            %Store the optimal choice for each state (by P);
            Inv0 = I0(ind0);
            Inv1 = I1(ind1);

            inv0_temp =  Inv0(:,:,:,ps);            
            Y0 = inv0_temp(:);
            estim_policy0{t,ps} = compute_policy_function(Y0, State_var);

            inv1_temp =  Inv1(:,:,:,ps);                        
            Y1 = inv1_temp(:);            
            estim_policy1{t,ps} = compute_policy_function(Y1, State_var);

        end

    else %Before Last period;

        grid_inv0 = inv0{t};
        grid_inv1 = inv1{t};


        [I0, hc0, H0, P_types0] =ndgrid(grid_inv0, hk{t}, grid_p, [0,1] );
        [I1, hc1, H1, P_types1] =ndgrid(grid_inv1, hk{t}, grid_p, [0,1] );

        %Distribution coming from the Solve_Forward;
        set_kids2 = fp_parfor.distrib_skills(:,t+1);
        set_PS2   = fp_parfor.distrib_PS(:,t);

        %Iteration
        set_kids = Hc_tp1_guess{t};
        set_kids1 = set_kids(:);

        %Random graph generation;
        [H1_tp1, H2_tp1] = ndgrid(set_kids1, set_kids2);

        set_PS_0 = zeros(length(set_kids1),1);
        set_PS_1 = ones(length(set_kids1),1);

        [PS1_0, PS2_0] = ndgrid(set_PS_0, set_PS2);
        [PS1_1, PS2_1] = ndgrid(set_PS_1, set_PS2);

        m_p0 = mean_peers(H1_tp1,H2_tp1,PS1_0,PS2_0,gamma,fp,set_kids2,fp_parfor,t);
        H_tp1_0 = m_p0.H_prime ;

        m_p1 = mean_peers(H1_tp1,H2_tp1,PS1_1,PS2_1,gamma,fp,set_kids2,fp_parfor,t);
        H_tp1_1 = m_p1.H_prime ;


        %Evolution of peer qualily by parenting style;
        H_prime_0 = reshape(nanmean(log(H_tp1_0),2),[size(hc0,1),size(hc0,2),size(hc0,3),size(hc0,4)]);       
        H_prime_1 = reshape(nanmean(log(H_tp1_1),2),[size(hc1,1),size(hc1,2),size(hc1,3),size(hc1,4)]);   
                                     
        %Evolution of a child's skills by parenting style;
        log_Skills_tp1_0 = log(technology(t,alpha, hc0, I0, 0, H0, TFP) ) ;
        log_Skills_tp1_1 = log(technology(t,alpha, hc1, I1, 1, H1, TFP) ) ;  

        %Calculate Altruism part of the parents' utility (based on equation 7);

        %Calculate Altruism based on Expected Value of Child's Utility (P=0);        
        temp_utility_child0 = nanmean( m_p0.U_links1,3);
        utility_child0 = nansum( temp_utility_child0,2);
        Altruism_0 = reshape( utility_child0 ,[size(hc0,1),size(hc0,2),size(hc0,3),size(hc0,4)]);


        %Calculate Altruism based on Expected Value of Child's Utility (P=1);        
        temp_utility_child1 = nanmean( m_p1.U_links1,3);
        utility_child1 = nansum( temp_utility_child1,2);
        Altruism_1 = reshape( utility_child1 ,[size(hc1,1),size(hc1,2),size(hc1,3),size(hc0,4)]);


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Solving the parent's problem;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;


        for ps = 1:1:2

            %Continuation value by parenting style (given parent type);
            Cont_value0(:,:,:,ps) = Interp_function(estim_value{t+1,ps},log_Skills_tp1_0(:,:,:,ps),H_prime_0(:,:,:,ps));
            Cont_value1(:,:,:,ps) = Interp_function(estim_value{t+1,ps},log_Skills_tp1_1(:,:,:,ps),H_prime_1(:,:,:,ps));

            %Bellman Equation (Equation 2) by Parenting Style and Type of parents (ps);     
            U0(:,:,:,ps) = fp.lambda.*delta(1).*log(hc0(:,:,:,ps)).*(1 - delta(3).*0) + (1-fp.lambda).*Altruism_0(:,:,:,ps) + 1.*log(1-I0(:,:,:,ps))  -(param_preference_PS_moved.*(P_types0(:,:,:,ps)==1)+delta(2).*(P_types0(:,:,:,ps)==0)).*0  + beta.*Cont_value0(:,:,:,ps);
            U1(:,:,:,ps) = fp.lambda.*delta(1).*log(hc1(:,:,:,ps)).*(1 - delta(3).*1) + (1-fp.lambda).*Altruism_1(:,:,:,ps) + 1.*log(1-I1(:,:,:,ps))  -(param_preference_PS_moved.*(P_types1(:,:,:,ps)==1)+delta(2).*(P_types1(:,:,:,ps)==0)).*1  + beta.*Cont_value1(:,:,:,ps);

        end


        [Value0,ind0] = max(U0,[],1);
        [Value1,ind1] = max(U1,[],1);

        temp_value0 = repmat( Value0 , [size(I0,1), 1, 1 , 1] );
        temp_value1 = repmat( Value1 , [size(I1,1), 1, 1 , 1] );


        %Optimal behavior about parenting style;
        I_guess0 = repmat(I0(ind0), [size(I0,1), 1, 1 , 1] );
        I_guess1 = repmat(I1(ind1), [size(I1,1), 1, 1 , 1] );

        Prob_PS = (1./(1+1./exp(temp_value1 -temp_value0)));

        %Updating the guess for next-period peers skills according to optimal behavior of parents;
        Hc_tp1_guess{t} = (1-Prob_PS).*technology(t,alpha, hc1, I_guess0, 0, H1, TFP) + Prob_PS.*technology(t,alpha, hc1, I_guess1, 1, H1, TFP);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Computing Value e Policy Functions;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;


        for ps = 1:1:2 %Loop over types of parents (moved or receiving);

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
            %Value Function;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
            Value0_actual(:,:,:,ps,t) = Value0(:,:,:,ps);
            Value1_actual(:,:,:,ps,t) = Value1(:,:,:,ps);


            %Collect State Variables;
            hc_s = hc0(1,:,:, ps);
            H_s  = H1(1,:,:, ps);
            State_var = [ hc_s(:), H_s(:) ];

            %Value for P=0;
            value0 = Value0(:,:,:,ps);           
            Y0 = value0(:);
            estim_value0{t, ps} = compute_value_function(Y0, State_var);    

            %Value for P=1;
            value1 = Value1(:,:,:,ps);
            Y1 = value1(:);
            estim_value1{t, ps} = compute_value_function(Y1, State_var);    

             %Expected Value over preference shocks;
            Y = log( exp(Y0) + exp(Y1));
            estim_value{t,ps}  = compute_value_function(Y, State_var);


            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Policy Function (Investments)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            %Store the optimal choice for each state (by P);
            Inv0 = I0(ind0);
            Inv1 = I1(ind1);

            inv0_temp =  Inv0(:,:,:,ps);            
            Y0 = inv0_temp(:);
            estim_policy0{t,ps} = compute_policy_function(Y0, State_var);

            inv1_temp =  Inv1(:,:,:,ps);                        
            Y1 = inv1_temp(:);            
            estim_policy1{t,ps} = compute_policy_function(Y1, State_var);

        end

    end  %End "if last period or not";

end %End loop for "t";

sb.Hc_tp1_guess_prime = Hc_tp1_guess;

sb.estim_policy0 = estim_policy0;
sb.estim_policy1 = estim_policy1;
sb.estim_value0  = estim_value0;
sb.estim_value1  = estim_value1;
sb.Value0_actual = Value0_actual;
sb.Value1_actual = Value1_actual;
